/*
    Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 
    Derivation ofthe xfer class for Metropolis Light Transport.
*/

#ifndef __mlt_xfer__
#define __mlt_xfer__

#include "xfer.h"
#include <iostream>
#include "boost/foreach.hpp"
#include "boost/lambda/lambda.hpp"
#include "boost/lexical_cast.hpp"
#include "threadlocal.h"
#include "emission.h"
#include <boost/math/special_functions/fpclassify.hpp>
#include <list>

using boost::lexical_cast;

namespace mcrx {
  template <typename> class vertex;
  template <typename> class ray_path;
  template<typename> class transition_return_type;
  template<typename T>
  std::ostream& operator<<(std::ostream& os, const mcrx::ray_path<T>& path);
  template<typename T> 
  std::ostream& compact(std::ostream& os, const blitz::Array<T,1>& a);
  template <typename, typename> class mlt_xfer;
}


template <typename T_lambda>
class mcrx::vertex {
public:
  /// The vertex location.
  vec3d pos_;
  /// The post-edge vertex factor (the opacity)
  T_lambda e_;
  /// The pre-edge vertex factor (the scattering function)
  T_lambda s_;
  /// The vertex sampling probability, as the ratio of p/s.
  T_float p_;

  vertex(const vec3d& pos, const T_lambda& e, const T_lambda& s, T_float p) :
    pos_(pos), e_(e), s_(s), p_(p) {};
  /// Constructor inits the vertex to be NaNs
  vertex(int n_lambda) :
    pos_(blitz::quiet_NaN(T_float())), e_(n_lambda), s_(n_lambda), 
    p_(blitz::quiet_NaN(T_float())) {
    e_=blitz::quiet_NaN(T_float()); s_=e_;};
  vertex(const vertex& v) :
    pos_(v.pos_), e_(independent_copy(v.e_)),
    s_(independent_copy(v.s_)), p_(v.p_) {};
};


template <typename T_lambda>
class mcrx::ray_path {
public:
  typedef vertex<T_lambda> T_vertex;

  /// The vertex quantities.
  std::list<T_vertex> v_;
  /// The edge factors (exp(-tau)/r^2)
  std::list<T_lambda> g_;

public:
  ray_path() {};
  ray_path(const ray_path& p) :
    // we can use default copy constructor for vertices since they
    // have a copy constructor defined
    v_(p.v_) {
    BOOST_FOREACH(const T_lambda& g, p.g_) {
      g_.push_back(make_thread_local_copy(g)); }};

  /// Returns the number of vertices in the path.
  int length() const {return v_.size(); };

  /** Checks whether the path is empty (i.e., only contains one
      vertex, which in that case must be a virtual vertex). */
  bool empty() const {
    assert (!v_.empty()); return (++v_.begin()) == v_.end(); };
  
  /** Checks the path direction by looking at the virtual vertices at
      the ends. 0 means it's a complete path, +1 means it's a source
      path, -1 means it's a camera path. */
  int direction() const {
    assert(!v_.empty());
    return 
      isnan(v_.front().e_(0)) ? 1 : 0 + 
      isnan(v_.back().s_(0)) ? -1 : 0; };


  /** Creates an empty ray which only contains the virtual
      endpoint. The direction argument indicates whether it's an empty
      source ray or an empty camera ray. */
  static ray_path create_empty(bool forward, int n_lambda) {
    ray_path p; 
    T_vertex v(vec3d(blitz::quiet_NaN(T_float())),
	       T_lambda(n_lambda), T_lambda(n_lambda), 1.0);
    (forward ? v.e_ : v.s_ ) = blitz::quiet_NaN(T_float());
    (forward ? v.s_ : v.e_ ) = 1.0;
    p.v_.push_back(v);
    T_lambda g(n_lambda);
    threadLocal_warn(g, true);
    g=1.0;
    p.g_.push_back(g);
    return p;
  }

  /**  Split the path at edge i, 0<=i<l. The edge i is removed, all
       vertices before (0-i) are still in the object, all vertices
       after (i+1 - n-1) are returned. If i<0, the entire path is
       copied to the return path. If i>=path length, the path is
       unchanged and an empty path is returned. */
  ray_path split(int i);
};


template<typename T_lambda>
mcrx::ray_path<T_lambda> mcrx::ray_path<T_lambda>::split(int i) 
{
  if(i<0) {
    // Move entire path to the return path.
    ray_path e;
    e.v_.splice(e.v_.begin(), v_, v_.begin(), v_.end());
    e.g_.splice(e.g_.begin(), g_, g_.begin(), g_.end());
    return e;
  }

  typename std::list<T_vertex>::iterator vi=v_.begin();
  typename std::list<T_lambda>::iterator gi=g_.begin();

  // Find split point by iterating from the front
  vi++;
  while(i--) {
    vi++;
    gi++;

    if(vi==v_.end()) {
      // i>=length, do nothing
      return ray_path();
    }
  }
  // XXX Must erase the appropriate s/e numbers for the vertex end
  // points that no longer make sense

  // vi is now pointing to the first vertex in the tail
  // gi is pointing to the edge to be erased
  gi=g_.erase(gi);
  // gi is now pointing to the first edge in the tail

  ray_path e;
  // splice the tails onto the new path
  e.v_.splice(e.v_.begin(), v_, vi, v_.end());
  if(gi!=g_.end())
    e.g_.splice(e.g_.begin(), g_, gi, g_.end());

  return e;
}

template<typename T> 
std::ostream& mcrx::compact(std::ostream& os,
			    const blitz::Array<T,1>& a)
{
  os << "(";
  for(int i=0;i<a.size();++i) {
    os << a(i);
    if(i<a.size()-1)
      os << ", ";
  }
  os << ")";
  return os;
}

template<typename T_lambda>
std::ostream& mcrx::operator<<(std::ostream& os, 
			       const mcrx::ray_path<T_lambda>& path)
{
  os << "\n[[ ";
  typename std::list<mcrx::vertex<T_lambda> >::const_iterator vi=path.v_.begin();
  typename std::list<T_lambda>::const_iterator gi=path.g_.begin();  
  int i=0;
  while(vi!=path.v_.end()) {
    os << i++ << " @ " << vi->pos_ << '\n';
    compact(os << "\te: ", vi->e_) << '\n';
    compact(os << "\ts: ", vi->s_) << "\tp: " << vi->p_ ;
    vi++;
    if(gi!=path.g_.end()) {
      compact(os << "\n\t\tg: ", *gi);
      gi++;
    }
    os << '\n';
  }
  os << "\t]]\n";
  return os;
}


template<typename T_lambda>
class mcrx::transition_return_type {
public:
  T_float p_ratio;
  T_float new_common;
  T_float old_common;
  T_lambda i_new;
  T_lambda i_old;
  boost::shared_ptr<ray_path<T_lambda> > path;
  std::string desc;

  transition_return_type(const boost::shared_ptr<ray_path<T_lambda> >& p, 
			 T_float nc, T_float oc,
			 const T_lambda& inew, 
			 const T_lambda& iold, T_float pr,
			 const std::string& d="") :
    path(p), new_common(nc), old_common(oc),
    i_new(make_thread_local_copy(inew)), 
    i_old(make_thread_local_copy(iold)), p_ratio(pr), desc(d) {};
  transition_return_type(const transition_return_type& rhs) :
    path(rhs.path), new_common(rhs.new_common), old_common(rhs.old_common),
    i_new(make_thread_local_copy(rhs.i_new)),
    i_old(make_thread_local_copy(rhs.i_old)), p_ratio(rhs.p_ratio),
    desc(rhs.desc) {};
};

template<typename T_iterator> void add_iter(T_iterator& i, int n) { while(n--)++i; };
template<typename T_iterator> void sub_iter(T_iterator& i, int n) { while(n--)--i; };


/** Responsible for performing the Monte-Carlo ray tracing.  The xfer
    class contains the machinery necessary to stochastically create
    rays, propagate these through the grid, scattering and/or
    absorbing the rays and finally detecting the rays leaving the grid
    volume. For a description of the algorithm used, see my thesis at
    http://sunrise.familjenjonsson.org/thesis . */
template <typename dust_model_type, typename grid_type>
class mcrx::mlt_xfer : public xfer<dust_model_type, grid_type> {
public:
  typedef xfer<dust_model_type, grid_type> T_base;

  typedef typename T_base::T_dust_model T_dust_model;
  typedef typename T_base::T_grid T_grid;
  typedef typename T_base::T_rng_policy T_rng_policy;
  typedef mcrx::T_float T_float;
  typedef typename T_base::T_scatterer T_scatterer;
  typedef typename T_base::T_lambda T_lambda;
  typedef typename T_base::T_biaser T_biaser;
  typedef typename T_base::T_ray T_ray;
  typedef typename T_base::T_emission T_emission;
  typedef typename T_base::T_emergence T_emergence;
  typedef typename T_base::T_cell T_cell;
  typedef typename T_base::T_cell_tracker T_cell_tracker;
  typedef typename T_base::T_propagation_history T_propagation_history;

  // typedefs for the MLT
  typedef ray_path<T_lambda> T_ray_path;
  typedef typename T_ray_path::T_vertex T_vertex;
  typedef typename std::list<T_vertex>::iterator T_vli;
  typedef typename std::list<T_lambda>::iterator T_lli;
  typedef typename std::list<T_vertex>::const_iterator T_vlci;
  typedef typename std::list<T_lambda>::const_iterator T_llci;

  int max_n_track_;
private:
  T_ray_path current_path_; ///<

  /// The accumulated probability of the current path (only used for the seed paths)
  T_lambda camera_intensity_;
  //std::vector<T_ray_path> seed_paths_;
  std::vector<T_float> seed_weights_;
  std::vector<std::pair<int,int> > seed_lengths_;
  //std::vector<T_float> seed_probs_;
  //std::vector<T_lambda> seed_intensities_;
  /// The emission distribution object
  boost::shared_ptr<const angular_distribution<T_lambda> > emdistr_;
  /// The camera distribution object
  boost::shared_ptr<const angular_distribution<T_lambda> > camdistr_;
  /// The dust (scatterer) distribution object
  boost::shared_ptr<const angular_distribution<T_lambda> > dustdistr_;

  T_float perturbation_scale_, image_perturbation_scale_, last_extra_weight_;
  T_float max_perturbation_scale_;
  T_float mutation_weight_, perturbation_weight_;

  void consolidate_arrays(array_1&);

public:
  mlt_xfer (/// The grid in which ray tracing is to be done
	    T_grid&,
	    /// The distribution of emission
	    const T_emission&,
	    /// The emergence object, with cameras detecting emerging rays.
	    T_emergence&,
	    /// The dust_model object
	    const T_dust_model&,
	    /// Seed for the random number generator.
	    T_rng_policy& rng,
	    /** The biaser to use. */
	    const T_biaser& b,
	    int minscat=0,
	    int maxscat=blitz::huge(int()));
  mlt_xfer (const mlt_xfer& rhs);
  ~mlt_xfer () {}

  T_lambda path_intensity(const T_ray_path&,  T_float& p_common, std::vector<T_float>& p_denom,
			  int start_vertex=0, int end_vertex=blitz::huge(int())) const;

  T_lambda path_intensity(const T_ray_path&,  
			  int start_vertex=0,
			  int end_vertex=blitz::huge(int())) const;
  T_float path_probability(const T_ray_path& p,
			   int start_vertex, int end_vertex) const;
  T_float total_path_probability(const T_ray_path& p,
				 int start_vertex, int end_vertex) const;
				 
  bool extend_path(T_ray_path& p, int n, int n_rr_start=blitz::huge(int()));
  void update_path_factors(T_ray_path& p, int start, int end);
  void merge_paths(T_ray_path& start_path, T_ray_path& end_path);
  int generate_bidir_paths(bool, bool,
			   int ns_max=blitz::huge(int()), 
			   int nc_max=blitz::huge(int()));

  T_ray_path recreate_bidir_path(std::pair<int,int> ns);

  T_ray_path make_path(const std::vector<vec3d>& pos);
  void pathtest();

  T_float sample_gaussian();
  T_float poisson(int i, T_float lambda) const;
  int sample_poisson(T_float) const;
  vec3d perturb_point(const vec3d& point, T_float scale=-1);
  T_float perturbation_scale(const vec3d& x, bool) const;
  T_float perturbation_probability(const vec3d& x0, const vec3d& x, 
				   T_float scale) const;
  transition_return_type<T_lambda> perturb_scattering_point(const T_ray_path& path);
  transition_return_type<T_lambda> bidir_mutation(const T_ray_path& path);

  /** Returns the deletion weight for the subpath between start_vertex
      and end_vertex. */
  T_float deletion_weight(const T_ray_path& p, 
			  int start_vertex, int end_vertex) const {
    const int n=end_vertex-start_vertex;
    assert(n>0);
    // pd1 only depends on the length
    const T_float pd1 = (n==1) ? 0.25 : ( (n==2) ? 0.5 : pow(2.0,-n) );
    // pd2 is the part of the acceptance probabilty that we can
    // calculate without knowing the proposal path. This is the total
    // path probability of the selected subpath/the path intensity of
    // the selected subpath
    const T_float pd2 = 
      total_path_probability(p, start_vertex, end_vertex);
    // we sum the intensity direct, which makes the compiler better
    // optimize out the allocation of the array completely.
    const T_float f_subpath=sum(path_intensity(p, start_vertex, end_vertex));
    //DEBUG(1,printf("\t\tDeletion weight: %5.3g\t%5.3g\t%5.3g\n",pd1,pd2,f_subpath););
    //const T_float pd2=1; const T_float f_subpath=1; // try without this
    const T_float w = (pd2>0) ? pd1*pd2/f_subpath : 0;
    return w; };


  /** Returns the cumulative deletion weights of all subpaths. */
  std::vector<T_float> total_deletion_weight(const T_ray_path& p) const {
    const int len=p.length();
    std::vector<T_float> dw;

    // loop over start of deleted subpath (last vertex to keep in s-part)
    //DEBUG(1,std::cout << "Subpath deletion weights:\n";);
    for(int s=1; s<=len-3; ++s)
      // loop over end of deleted subpath
      for(int t=s+1; t<=len-2; ++t) {
	const T_float w = deletion_weight(p, s, t);
	dw.push_back(w + ((dw.size()==0) ? 0: dw.back() ));
	//DEBUG(1,printf("\t%d\t%d\t%5.3g\n",s,t,w););
      }
    return dw;
  };

  /** Returns the weight of adding n vertices during a mutation, based
      on the fact that n_d were deleted. len is the original length of
      the path. We do a special case in case no vertex was deleted so
      that we then always add at least one. If we're not including the
      direct path, we also always add at least one if the deleted
      points would make a path without scattering vertices. */
  T_float addition_weight(int len, int n_d, int n, 
			  bool include_direct) const {
    const int dn = n; //n - n_d;
    bool must_add = (n_d==0);
    if(!include_direct) must_add |= (len-n_d==4);

    return (dn==0) ? (must_add ? 0 : 0.25) : ( (dn==1) ? 0.15 : 0.2*pow(2.0,-dn) );
  };

  /** Returns the total addition weight summed over the possible n>=0,
      when n_d were deleted. len is the original length of the path. */
  T_float total_addition_weight(int len, int n_d, bool include_direct) const {
    // the n>=n_d is an infinite sum that can be done
    // analytically. sum( 0.5^i )=2, but we are only summing for i>=2,
    // so that sum is 0.5. The total is 0.75.
    bool must_add = (n_d==0);
    if(!include_direct) must_add |= (len-n_d==4);

    T_float sum = (must_add ? 0 : 0.25) + 0.15 + 0.2*0.5;
    // the n<n_d we must do explicitly
    //for(int n=n_d-1; n>=0; --n)
    //sum+=addition_weight(n_d, n);
    return sum;
  };
    
  /** Returns the probability weight of perturbing a vertex. */
  T_float perturbation_weight(const T_ray_path& path) const {
    // can only mutate if we have non-endpoints
    return path.length()>4 ? 
      perturbation_weight_ :
      0.0;
  };

  /** Probability of perturbing vertex v. */
  T_float p_perturb_vertex(const T_ray_path& p, int v) const {
    const int n = p.length()-4;
    assert(n>0);
    const T_float extra_weight=last_extra_weight_*n;
    const T_float prob = ((v==n+1) ? 1.+extra_weight : 1.)/(n+extra_weight);
    return prob;
  };

  /** Returns the probability weight of doing a bidir mutation. */
  T_float mutation_weight(const T_ray_path& p) const {
    return mutation_weight_; };

  /** Returns the total probability weight of mutation possibilities
      for the specified path. */
  T_float total_weight(const T_ray_path& p) const {
    return mutation_weight(p) + perturbation_weight(p);
  };

  /** Return probability of doing a perturbation. */
  T_float p_perturbation(const T_ray_path& p) const {
    return perturbation_weight(p)/total_weight(p);};

  /** Return probability of a budirectional mutation. */
  T_float p_mutation(const T_ray_path& p) const {
    return mutation_weight(p)/total_weight(p);};

  transition_return_type<T_lambda> proposal(const T_ray_path& path);

  int metropolis(const T_ray_path& start_path, int n, T_float w0,
		 bool write_debug);
  T_float metropolis_driver(int nseed, int nmetropolis, int nmutations);
  void bidirectional_shoot(int ns_max, int nc_max);
  const T_ray_path& current_path() const {return current_path_;};
};


/** Recreates a path with specific n,s generated by
    generate_bidir_paths, assuming the RNG state is identical before
    the call. */
template <typename dust_model_type, typename grid_type>
typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_ray_path
mcrx::mlt_xfer<dust_model_type, grid_type>::
recreate_bidir_path(std::pair<int,int> ns_len)
{
  const int n_lambda = this->norm_.size();
  T_ray_path source_path(T_ray_path::create_empty(true, n_lambda));
  T_ray_path camera_path(T_ray_path::create_empty(false, n_lambda));

  extend_path(source_path, blitz::huge(int()), 2);
  extend_path(camera_path, blitz::huge(int()), 2);
  const int ns=source_path.length();
  const int nc=camera_path.length();
  const int t = ns_len.first - ns_len.second;
  T_ray_path s_subpath(source_path); s_subpath.split(ns_len.second+1);
  T_ray_path c_subpath(T_ray_path(camera_path).split(nc-t-3));
  merge_paths(s_subpath, c_subpath);
  return s_subpath;
}

/** Generates a number of paths by bidirectional sampling. Two paths
    are generated and then joined in all possible ways, adding them to
    the camera. */
template <typename dust_model_type, typename grid_type>
int
mcrx::mlt_xfer<dust_model_type, grid_type>::
generate_bidir_paths(bool add_to_camera, bool ignore_direct,
		     int ns_max, int nc_max)
{
  assert(ns_max>=1);
  assert(nc_max>=1);

  const int n_lambda = this->norm_.size();
  DEBUG(1,std::cout << "\nGenerating bidirectional path"<<std::endl;);

  T_ray_path source_path(T_ray_path::create_empty(true, n_lambda));
  T_ray_path camera_path(T_ray_path::create_empty(false, n_lambda));

  // Extend both paths until they leave the volume
  //const int ns_max=blitz::huge(int());
  //const int nc_max=blitz::huge(int()); //1; // emulate normal Sunrise
  // note, if neither ns_max or nc_max == inf, the results will be
  // biased because we'll be truncating the possible path space
  extend_path(source_path, ns_max, 2);
  extend_path(camera_path, nc_max, 2);
  DEBUG(1,std::cout << "Source path length " << source_path.length() << ", camera path length " << camera_path.length() << std::endl;);

  DEBUG(2,std::cout << "Source path: "<< source_path << std::endl;);
  DEBUG(2,std::cout << "Camera path: "<< camera_path << std::endl;);

  // Now we generate complete paths by merging all subpaths of the
  // source path with all subpaths of the camera path.
  const int ns=source_path.length();
  const int nc=camera_path.length();
  int n_total=0;

  // first loop is over the number of *added* vertices, which can be a
  // maximum of ns+nc-4
  for(int n=(ignore_direct?1:0); n<=ns+nc-4; ++n) {
    // Now we need to loop over the different ways to add n vertices.
    // There are up to n+1 ways, but this is limited by the maximum
    // length of the source path, so the realizable ways are ns-1.

    // The loop variable is number of vertices taken from the camera
    // path. There are a maximum of ns-2 vertices to take from the
    // source path, and we can take at most nc-2 from the camera path,
    // so the minimum we need from the source path is n-nc+2.
    const int s_min = std::max(0, n-nc+2);
    const int s_max = std::min(n, ns-2);
    DEBUG(1,std::cout << "\n\nProcessing paths with " << n << " added vertices\n";);
    DEBUG(2,std::cout << "Taking " << s_min <<"-" << s_max << " from source path\n";);

    for(int s=s_min; s<=s_max; ++s) {

      const int t = n-s;

      // subpath lengths are s+2,t+2
      DEBUG(1,printf("\nAdded vertices: %d,%d\n",s,t););

      T_ray_path s_subpath(source_path); s_subpath.split(s+1);
      T_ray_path c_subpath(T_ray_path(camera_path).split(nc-t-3));

      //DEBUG(2,std::cout << "Source subpath: " << s_subpath << std::endl;);
      //DEBUG(2,std::cout << "Camera subpath: " << c_subpath << std::endl;);
      
      merge_paths(s_subpath, c_subpath);
      // merged path has length n+4 = s+t+4

      DEBUG(2,std::cout << "\nMerged path: " << s_subpath;);

      // common is the factor common to both the intensity and the
      // probabilities. any ratio of probabilities and intensities is
      // thus independent of it.
      T_float common;
      std::vector<T_float> p_edge;
      T_lambda cur_i(path_intensity(s_subpath, common, p_edge));

      DEBUG(1,compact(std::cout << "\nPath intensity ", cur_i) << std::endl;);
      DEBUG(1,std::cout << "Common factor " << common << std::endl;);

      T_float w;
      T_float prob;

      const bool use_balance_heuristic=true;
      if(use_balance_heuristic) {
	// the balance heuristic is to calculate the contribution as
	// f/(sum p), where the sum is the total probability of drawing
	// this sample using all the techniques. The weighting factor is
	// p/(sum p). We thus need to loop again over the different
	// possible ways we *could have* generated this sample. Note
	// that this does not depend on the actual number of ns and nc,
	// because if it wasn't realized that's already taken care of by
	// the probability. It does however depend on the ns_max and
	// nc_max.

	DEBUG(2,std::cout << "Evaluating subpath probabilities:\n";);
	// ss is the number of vertices that "could have been" taken from the source path
	DEBUG(2,printf("\tC-edge\tFactor\tProb\n"););
      
	T_float ptot=0;
	const int ss_min = std::max(0, n-nc_max+1);
	const int ss_max = std::min(n, ns_max-1);
	for(int ss=ss_min; ss<=ss_max; ++ss) {
	  const int tt = n-ss;
	  //const T_float pp1 = path_probability(s_subpath, 1, 1+ss);
	  //const T_float pp2 = path_probability(s_subpath, n+2, n+2-tt);
	  ptot += p_edge[ss];
	  if(ss==s)
	    prob = p_edge[ss];

	  DEBUG(2,printf("%s\t%d-%d\t%5.3g\t%5.3g\n",(ss==s)?"*":"",1+ss ,2+ss, p_edge[ss], common*p_edge[ss]));
	}

	w=prob/ptot;
      }
      else {
	// simplest (and dumbest) weighting is to weight all different
	// s combinations by 1/n+1
	prob = p_edge[s];
	//path_probability(s_subpath, 1, 1+s) *
	//path_probability(s_subpath, n+3, n+3-t);
	w=1./(s_max-s_min+1);
      }

      DEBUG(1,printf("Weight: %g\n", w));

      if(add_to_camera) {
	if(!(w/prob < blitz::huge(T_float()))) {
	  DEBUG(1,std::cout << "\tProbability underflow - skipping path\n";);
	  continue;
	}
	cur_i *= w/prob;
	DEBUG(1,compact(std::cout << "Image contribution ", cur_i) << std::endl;);
	if(any(cur_i)>1e25) {
	  std::cout << "LARGE IMAGE CONTRIBUTION " << cur_i << std::endl;
	  std::cout << "Path:" << s_subpath << std::endl;
	}
	if( n>=this->n_scat_min_ && n<=this->n_scat_max_) {
	  blitz::TinyVector<int,2> px =
	    (*this->current_emergence_->begin())->project((++(++s_subpath.v_.rbegin()))->pos_);
	  if(px[0]>=0) {
	    // put contribution in correct camera if we are disaggregating
	    if(max_n_track_<0)
	      (*this->current_emergence_->begin())->add(px, 1., cur_i);
	    else {
	      const int cam= (n<=max_n_track_)? n*(n+1)/2+s : 
		(max_n_track_+1)*(max_n_track_+2)/2;
	      assert(cam<this->current_emergence_->n_cameras());
	      (*(this->current_emergence_->begin()+cam))->add(px, 1., cur_i);
	    }
	  }
	}
      }
      else {
	// in this case we are generating seed paths for metropolis
	// runs, so save path weight and configuration. note that the
	// probability to save is the *total* probability of
	// generating the path, so prob/w.
	T_float cur_weight = (!(common*prob/w>0) || all(cur_i==0)) ? 
	  0.0 : sum(cur_i)*(w/prob);
	assert(cur_weight<blitz::huge(T_float()));
	seed_weights_.push_back(cur_weight + 
				((seed_weights_.empty()) ? 
				 0.0 : 
				 seed_weights_.back()));
	seed_lengths_.push_back(std::make_pair(n,s));
	assert(any(cur_i>0)||(cur_weight==0));
      }
      
      n_total++;
    }
  }
  DEBUG(1,std::cout << "\nAll subpaths processed.\n\n";);
  return n_total;
}


/** Calculates the total probability of generating the specified
    subpath between start_vertex and end_vertex by summing over all
    (s,t) combinations. (We *could* do this more effiiently by using
    the p_factor vector, but that will only give the probablity up to
    a common factor, because the p_common is appropriate for the
    entire path.) */
template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
total_path_probability(const T_ray_path& p,
		       int start_vertex, int end_vertex) const
{
  // number of (internal) vertices on subpath
  const int n=end_vertex-start_vertex-1;
  assert(n>=0);

  // when evaluating the total probability, we need to sum over all
  // combinations of adding the n vertices. There are n+1 such
  // combinations.
  T_float ptot=0;
  for(int s=0; s<=n; ++s) {
    const int t = n-s;
    const T_float pp1 = path_probability(p, start_vertex, start_vertex+s);
    const T_float pp2 = path_probability(p, end_vertex, end_vertex-t);
    ptot += pp1*pp2;
    
    //DEBUG(1,printf("\t\t%s\t%d-%d\t%d-%d\t%5.3g\t%5.3g\t%5.3g\n","",start_vertex ,start_vertex+s, end_vertex,end_vertex-t, pp1,pp2,pp1*pp2));
  }
  //DEBUG(1,printf("\t\tTotal %d-%d path probability: %5.3g\n",start_vertex,end_vertex,ptot););
  assert(ptot<blitz::huge(T_float()));

  return ptot;
}



/** Extends a (possibly empty, i.e. only containing the virtual start
    vertex for either the source or the camera) path with n vertices
    by casting rays using the normal MC algorithm. Whether the path is
    extended from the front or back is determined by the path
    direction. If n vertices successfully were added, true is
    returned. If the ray exited at some point, false is returned and
    the path contains the vertices that were added.  n_rr_start
    determines after how many vertices russian roulette (based on the
    albedo at the vertex) will start. */
template <typename dust_model_type, typename grid_type>
bool
mcrx::mlt_xfer<dust_model_type, grid_type>::
extend_path(T_ray_path& p, int n, int n_rr_start)
{
  if (n==0) 
    return true;

  const int dir = p.direction();
  const int n_lambda = this->norm_.size();

  assert(dir!=0);
  assert(n>0);

  T_float prob;
  bool add_end = p.length()<=2;
  vec3d original_direction;
  T_float intensity=1;

  // there are 3 possibilities here: 1. the path is "empty", in which
  // case it only has a virtual vertex and the source/camera vertex
  // needs to be added. 2. the path has a source/camera vertex but
  // nothing more. In this case we need to redraw the *direction* out
  // from the source/camera but not the location. 3. The path ends
  // with a real scattering vertex, in which case we just do a
  // scattering
  if(add_end) {
    if(p.empty()) {
      /* If the path is "empty" it only has the virtual vertex and
	 edge. Then we need to add the source/camera vertex position and
	 e-factor.     */

      // This iterator points to the only vertex we have
      T_vli vi = (dir>0) ? (--p.v_.end()) : p.v_.begin();

      ray_distribution_pair<T_lambda> e = 
	dir>0 ? 
	this->emission_.emit(this->b_, this->rng(), prob) : 
	(*this->current_emergence_->begin())->emit(this->b_, this->rng(), prob);

      this->ray_ = *e.ray_pointer();
      
      T_vertex v(n_lambda);
      v.pos_ = this->ray_.position();
      v.e_ = 1;
      ((dir>0) ? vi->s_ : vi->e_) = e.normalization();
      // s-factor and probablitity are set after entering main loop
      if(dir>0) 
	p.v_.push_back(v); 
      else 
	p.v_.push_front(v); 
      n--;
    }
    else {
      // if the path has the source/camera vertex position, then we
      // just redraw the direction, which is a "scatter". (This
      // distinction is only meaningful for distributed emissions (and
      // non-pinhole cameras))

      // for this to work correctly on non-pointsource emissions, we
      // need to implement the scatter function for emitters.
      this->ray_.set_position((dir>0) ? p.v_.back().pos_ : p.v_.front().pos_);
      if(dir>0) {
	ray_distribution_pair<T_lambda> e = 
	  this->emission_.emit(this->b_, this->rng(), prob);
	this->ray_.set_direction(e.ray_pointer()->direction());
      }
      else {
	(*this->current_emergence_->begin())->scatter(this->b_, this->ray_, 
						      this->rng(), prob);
      }
    }

    // set up the (fake) original direction
    original_direction = vec3d(blitz::quiet_NaN(T_float()));
  }
  else {
    /* If the edge is not empty, we do the first scattering (by
       analogy with emission) here. */
    assert(!isnan(p.v_.back().e_(0)));

    // This iterator points to the vertex BEFORE the last
    T_vli vi = (dir>0) ? (--(--p.v_.end())) : (++p.v_.begin());

    // Initialize ray according to last edge
    this->ray_.set_position((dir>0) ? p.v_.back().pos_ : p.v_.front().pos_);
    vec3d dr= (this->ray_.position() - vi->pos_);
    T_float l=sqrt(dot(dr,dr));
    this->ray_.set_direction(dr/l);
    // and scatter it
    original_direction = this->ray_.direction();
    this->dust_model.get_scatterer(0).scatter(this->b_, this->ray_, 
					      this->rng(), prob);
  }

  // At this point we have a valid end vertex and the ray is set up to
  // start tracing in the current direction
  this->ray_.set_traversed_column(this->dust_model.zero_density() );
  this->current_cell = this->g.locate(this->ray_.position(), false);  

  // start looping over remaining vertices.
  while(n--) {

    // first calculate the angular distribution of the vertex we just
    // scattered
    const angular_distribution<T_lambda>& distr = 
      (add_end) ? 
      ( (dir>0) ?
	*emdistr_.get() :
	*camdistr_.get() ) :
      *dustdistr_.get();

    T_lambda d;
    assign(d, distr.distribution_function(this->ray_.direction(), 
					  original_direction));

    // if we used the scatterer, we also need to apply the albedo...
    T_float rr_factor=1;
    if(!add_end) {
      d *= this->dust_model.get_scatterer(0).albedo();

      // optionally, do a Russian roulette based on the albedo to see
      // if we should continue
      if(n_rr_start-- < 0) {
	const double P_rr=sum(this->dust_model.get_scatterer(0).albedo())/n_lambda;
	if(this->rng().rnd() < P_rr) {
	  // survived, decrease the probability (this will boost the
	  // contribution to make it unbiased)
	  rr_factor=P_rr;
	}
	else {
	  // died, kill it.
	  return false;
	}
      }
    }

    add_end=false;

    // Draw an optical depth and propagate
    this->scattering_reftau_ = - log (this->rng().rnd ());
    assert(this->scattering_reftau_>0);
      
    // Propagate ray until scatter or exit
    if(!this->current_cell) {
      this->propagate_external();
    }

    if(this->current_cell)
      // for some reason, if we don't supply an argument here, it
      // doesn't compile
      while( !(this->propagate<true, false, false, false>(blitz::huge(T_float()))) &&
	     this->current_cell);

    if(!this->current_cell)
      // we exited the volume (or we never hit the grid at all), which means path extension failed
      return false;

    assert(this->ray_.length()>0);
    this->scattering_tau_ =
      this->dust_model.total_abs_length(this->ray_.traversed_column()(blitz::tensor::j));

    // Add info. Need to set s- and p-factors of the last vertex, then
    // add the g-factor of the edge (exp-tau/r^2) and the position and
    // e-factor of new vertex (kappa)

    T_vertex v(n_lambda);
    v.pos_ = this->ray_.position();
    v.e_ = this->dust_model.total_abs_length(this->current_cell.data()->get_absorber().densities()(blitz::tensor::j));

    if(dir>0) {
      p.v_.back().s_ = d;
      // set the p-factor to the probability over the distribution reference
      p.v_.back().p_ = (prob>0) ? prob*rr_factor/this->b_.reference1(d) : 0;

      p.g_.push_back(make_thread_local_copy(exp(-this->scattering_tau_)/
					    (this->ray_.length()*this->ray_.length())));
      p.v_.push_back(v);
    }
    else {
      p.v_.front().s_ = d;
      // set the p-factor to the probability over the distribution reference
      p.v_.front().p_ = (prob>0) ? prob*rr_factor/this->b_.reference1(d) : 0;

      p.g_.push_front(make_thread_local_copy(exp(-this->scattering_tau_)/
					     (this->ray_.length()*this->ray_.length())));
      p.v_.push_front(v);
    }

    // finally, we need to do the scattering associated with the
    // current vertex (for the last vertex this is never used but
    // that's ok). We know this is a dust scattering.
    original_direction = this->ray_.direction();
    this->dust_model.get_scatterer(0).scatter(this->b_, this->ray_, 
					      this->rng(), prob);

    // That completes this iteration.
  }

  return true;
}


/* Merges two paths by splicing the end_path onto the end of the
   start_path. The path is updated with the factors for the added
   edge. */
template <typename dust_model_type, typename grid_type>
void
mcrx::mlt_xfer<dust_model_type, grid_type>::
merge_paths(T_ray_path& start_path, T_ray_path& end_path) 
{
  assert(start_path.direction()>0);
  assert(end_path.direction()<0);

  const int n = start_path.length();

  // First add an edge to the end of the start path. We copy the last edge.
  start_path.g_.push_back(independent_copy(start_path.g_.back()));
  
  // Now splice the end_path onto it
  start_path.v_.splice(start_path.v_.end(), end_path.v_,
		       end_path.v_.begin(), end_path.v_.end());
  start_path.g_.splice(start_path.g_.end(), end_path.g_,
		       end_path.g_.begin(), end_path.g_.end());

  update_path_factors(start_path, n-1, n);
}


/** Updates the vertex and edge factors when a path has been
    changed. (In some other way than with extend_path, since it
    calculates the factors directly.) It only operates on the vertices
    start-end (inclusive) and the edges between them. */
template <typename dust_model_type, typename grid_type>
void
mcrx::mlt_xfer<dust_model_type, grid_type>::
update_path_factors(T_ray_path& p, int start, int end)
{
  const int len=p.length();

  assert(start>=1);
  assert(end<=len-2);
  assert(start<end);

  // Get iterators. We keep 2 running vertex iterators so we can
  // calculate directions
  T_vli v1 = p.v_.begin();
  add_iter(v1, start-1);
  T_vli v0 = v1++;
  T_lli gi = p.g_.begin();
  add_iter(gi, start);

  // calculate previous edge direction
  vec3d original_direction;
  original_direction = v1->pos_ - v0->pos_;
  original_direction /= mag(original_direction);
  ++v0; ++v1;
  int i=start;

  while(i<end) {

    // set up ray to start of edge
    this->ray_.set_position(v0->pos_);
    const vec3d dr= (v1->pos_ - v0->pos_);
    const T_float l=mag(dr);
    this->ray_.set_direction(dr/l);
    this->ray_.set_traversed_column(this->dust_model.zero_density() );

    // update s-factor for vertex. 
    // pick correct distribution object
    const angular_distribution<T_lambda>& distr = 
      (i==1) ? 
      *emdistr_.get() :
      ( (i==len-2) ?
	*camdistr_.get() :
	*dustdistr_.get() );
    assert(i!=len-2); // if this is true we use incorrect direction for ray
    v0->s_ = distr.distribution_function(this->ray_.direction(), original_direction);
    if((i>1)&&(i<len-2))
      // for the dust nodes, we also need to apply the albedo
      v0->s_ *= this->dust_model.get_scatterer(0).albedo();

    // update p-factor. XXX we ASSUME here that the probability is the
    // distribution function. We don't have an independent way of
    // getting it here currently
    const T_float prob = 
      this->b_.reference1(distr.distribution_function(this->ray_.direction(), 
						     original_direction));
    v0->p_ = (prob>0) ? prob/this->b_.reference1(v0->s_) : 0.0;

    // integrate the edge to get edge factor
    this->current_cell = this->g.locate(this->ray_.position(), false);
    if(!this->current_cell)
      this->propagate_external (l);
    
    // propagate to next point in path. propagate will return true if
    // we hit the point (but ALSO set current_cell to zero)
    bool scattered;
    while(this->current_cell &&
	  !(scattered= this->propagate<false, false, false, false>(l)) ) {}

    // if we come out here and scattered is false, we left the
    // volume. otherwise we hit the point

    // Add the g-factor of the edge (exp-tau/r^2) and the e-factor of
    // new vertex (kappa, unless it's the camera vertex in which case it's 1)
    *gi = exp(-this->dust_model.total_abs_length(this->ray_.traversed_column()(blitz::tensor::j)))/(l*l);
    
    if (i<len-3) {
      if(this->current_cell)
	v1->e_ = this->dust_model.total_abs_length(this->current_cell.data()->get_absorber().densities()(blitz::tensor::j));
      else
	// if we don't have a cell we know that the kappa is zero.
	v1->e_ = 0.0;
    }
    else
      v1->e_ = 1.0;

    // that completes this edge. 
    original_direction = this->ray_.direction(); 
    ++v0; ++v1;
    ++gi;
    ++i;
  }
  // when we come out here, we still need to update the s- and p-factors for
  // the end vertex.
  vec3d end_dir= (v1->pos_ - v0->pos_);
  end_dir /= mag(end_dir);

  // pick the correct distribution. the camera distribution requires
  // special handling here because it takes the outgoing direction
  // just like the emitters, but we are traversing the path in a
  // forward manner.
  if (i==len-2) {
    const angular_distribution<T_lambda>& distr = *camdistr_.get();
    v0->s_ = distr.distribution_function(-original_direction, end_dir);
  }
  else {
    const angular_distribution<T_lambda>& distr = 
      (i==1) ? 
      *emdistr_.get() :
      *dustdistr_.get();
    v0->s_ = distr.distribution_function(end_dir, original_direction);
  }
  // XXX ASSUME probabilty is the distribution function
  const T_float prob = this->b_.reference1(v0->s_);

  if((i>1)&&(i<len-2))
    // if it's a scattering vertex, we also need to apply the albedo
    v0->s_ *= this->dust_model.get_scatterer(0).albedo();
  v0->p_ = (prob>0) ? prob/this->b_.reference1(v0->s_) : 0;
}

/** Calculates the intensity and probability associated with a
    path. Because both the intensity and the probability largely have
    the same factors in them, it is advantageous to factor out the
    common factor. Because the relevant quantity (generally) is the
    intensity/probability, in RFOT situations both can underflow even
    though the weight of the path is not small. Factoring out the
    common factor avoids this behavior. So, to get the actual path
    intensity, the returned array should be multiplied by the common
    factor. The probability is such that the probability of the path
    depends on the connecting subpath, so to get the probability of
    generating a path through bidirectional tracing, you take the
    common factor multiplied by the edge_factor for the connecting
    edge. */
template <typename dust_model_type, typename grid_type>
typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_lambda
mcrx::mlt_xfer<dust_model_type, grid_type>::
path_intensity(const T_ray_path& p, T_float& common, 
	       std::vector<T_float>& p_edge,
	       int start_vertex, int end_vertex) const
{
  // Since all the data are in the ray path now, this is a simple
  // function. The intensity includes both the vertex factors on the start
  // and end vertices. To get the whole path, specify 0,HUGE_INT.
  p_edge.clear();
  const int n_lambda = this->norm_.size();

  if(p.v_.empty()) {
    common=0;
    return initialize(this->norm_, 0.0);
  }

  T_vlci vi=p.v_.begin();
  // v1 points to the *next* vertex, for probability calc
  T_vlci v1=++p.v_.begin();
  T_llci gi=p.g_.begin();
  T_lambda f(initialize(this->norm_, 1.0));

  // for the probability, we need to factor out the luminosity which
  // is in the s-factor for vertex 0. not any more
  common=1;
  T_float p_common=1./this->b_.reference1(vi->s_);
  int i=0;
  T_float lum;

  while(gi!=p.g_.end()) {
    T_float temp_ref;
    if(i==end_vertex) {
      temp_ref = this->b_.reference1(vi->e_);
      //assert(temp_ref<blitz::huge(T_float()));
      f *= vi->e_/temp_ref;
      common *= temp_ref;
      //assert((common==0) || all(f<blitz::huge(T_float())));
      // exclude s-factor if it's the end of the entire path, which is NaN.
      if(v1!=p.v_.end()) {
	temp_ref = this->b_.reference1(vi->s_);
	//assert(temp_ref<blitz::huge(T_float()));
	f *= vi->s_/temp_ref;
	common *= temp_ref;
	//assert((common==0) || all(f<blitz::huge(T_float())));
      }
    }
    else if((i>=start_vertex)&&(i<end_vertex)) {
      T_lambda ss(vi->s_);
      T_lambda gg(*gi);
      T_lambda oldf(independent_copy(f));
      T_float oldcommon=common;
      temp_ref = this->b_.reference1(vi->s_)*this->b_.reference1(*gi);
      //assert(temp_ref<blitz::huge(T_float()));
      f *= vi->s_ * *gi/temp_ref;
      common *= temp_ref;
      //assert((common==0) || all(f<blitz::huge(T_float())));

      // exclude e-factor from vertex 0, since it's NaN.
      if(i>0) {
	temp_ref = this->b_.reference1(vi->e_);
	//assert(temp_ref<blitz::huge(T_float()));
	f *= vi->e_/temp_ref;
	common *= temp_ref;
	//assert((common==0) || all(f<blitz::huge(T_float())));
      }
    }

    if((i>0) && (v1!=p.v_.end())) {
      const T_float pp= vi->p_;
      p_common *= pp;
      const T_float a=this->b_.reference1( vi->s_ );
      const T_float c=this->b_.reference1( *gi );
      const T_float d=this->b_.reference1( v1->s_ );
      const T_float e= v1->p_;
      const T_float pp_edge = 1.0/(a*pp*c*d*e);
      //assert(pp_edge==pp_edge);
      p_edge.push_back(pp_edge);
    }

    vi++;
    if(v1!=p.v_.end())
      v1++;
    gi++;
    i++;
    assert(vi!=p.v_.end());
  }

  // the common probability factor needs to be be factored into each
  // of the edge probabilities.
  using boost::lambda::_1;
  transform(p_edge.begin(), p_edge.end(), p_edge.begin(), _1*p_common);

  // when we get out here we still have the last vertex (since there
  // is one more vertex than edge. But there's no info in it.
  assert(start_vertex<i);

  // if common==0, we probably have the result that f contains some
  // nan or inf. in that case, we set it to zero.  Also, if
  // common=nan, it has probably increased to inf and then been
  // multiplied by zero. in that case, it's also zero.
  if((common<blitz::tiny(T_float())) || (common!=common)) {
    f=0;
    common=0;
  }
  assert(all(f==f));

  return f;
}


template <typename dust_model_type, typename grid_type>
typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_lambda
mcrx::mlt_xfer<dust_model_type, grid_type>::
path_intensity(const T_ray_path& p,
	       int start_vertex, int end_vertex) const
{
  // Since all the data are in the ray path now, this is a simple
  // function. The intensity includes both the vertex factors on the start
  // and end vertices. To get the whole path, specify 0,HUGE_INT.

  const int n_lambda = this->norm_.size();

  T_vlci vi=p.v_.begin();
  // v1 points to the *next* vertex, for probability calc
  T_vlci v1=++p.v_.begin();
  T_llci gi=p.g_.begin();
  T_lambda f(initialize(this->norm_, 1.0));
  f=1;

  int i=0;
  T_float lum;

  while(gi!=p.g_.end()) {
    if(i==end_vertex) {
      f *= vi->e_;
      // exclude s-factor if it's the end of the entire path, which is NaN.
      if(v1!=p.v_.end())
	f *= vi->s_;
    }
    else if((i>=start_vertex)&&(i<end_vertex)) {
      f *= vi->s_ * *gi;
      // exclude e-factor from vertex 0, since it's NaN.
      if(i>0)
	f *= vi->e_;
    }

    vi++;
    if(v1!=p.v_.end())
      v1++;
    gi++;
    i++;
    assert(vi!=p.v_.end());
  }

  // when we get out here we still have the last vertex (since there
  // is one more vertex than edge. But there's no info in it.
  assert(start_vertex<i);
  assert(all(f==f));
  return f;
}


template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
path_probability(const T_ray_path& p,
		 int start_vertex, int end_vertex) const
{
  // Since all the data are in the ray path now, this is a simple
  // function. The path probability includes the s-factor*p-factor on
  // the start vertex and the e-factor on the end vertex. The
  // direction matters for the sampling probability, so end can be <
  // start (but either may not include the end points which doesn't make sense)

  // it's convenient to say that p=1 if start=end
  if(start_vertex==end_vertex) return 1;

  assert(start_vertex>=0);
  assert(end_vertex>=0);
  assert(start_vertex<=p.v_.size()-1);
  assert(end_vertex<=p.v_.size()-1);

  const bool forward = (start_vertex<end_vertex);
  T_vlci vi= p.v_.begin();
  T_llci gi= p.g_.begin();

  T_float prob=1;
  int i=0;

  while(gi!=p.g_.end()) {

    if(i==start_vertex) {
      const T_float sp = this->b_.reference1(vi->s_) * vi->p_;
      prob *= sp;
      // this edge is included for forward traversals
      if(forward) {
	const T_float g = this->b_.reference1(*gi);
	prob *= g;
      }
    }
    if(i==end_vertex) {
      const T_float e = this->b_.reference1(vi->e_);
      prob *= e;
      // this edge is included for backward traversals
      if(!forward) {
	const T_float g = this->b_.reference1(*gi);
	prob *= g;
      }
    }
    if((i-start_vertex)*(i-end_vertex)<0) {
      const T_float seg = this->b_.reference1(vi->s_*vi->e_* *gi);
      const T_float p = vi->p_;
      prob *= seg * p;
    }

    vi++;
    gi++;
    i++;
  }

  return prob;
}


inline mcrx::T_float gaussian(mcrx::T_float x)
{
  return sqrt(2/constants::pi)*exp(-0.5*x*x);
}


template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
sample_gaussian()
{
  T_float x;
  while(x=this->rng().rnd()*5,this->rng().rnd()>gaussian(x)){};
  return x;
}


template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
poisson(int i, T_float lambda) const
{
  T_float p = pow(lambda, i)*exp(-lambda);
  // compute factorial (inefficient but simple...)
  T_float fact=1;
  for (int j=2; j<=i; ++j)
    fact*=j;
  return p/fact;
}


/** Poisson sampling due to Knuth. */
template <typename dust_model_type, typename grid_type>
int
mcrx::mlt_xfer<dust_model_type, grid_type>::
sample_poisson(T_float lambda) const
{
  T_float L=exp(-lambda), p=1;
  int k=0;
  do {
    k++;
    p*=this->rng_.rnd(); 
  }
  while (p>L);
  return k-1;
}


/** Perturbs a point by an exponential distribution distance in a
    random direcotion. If scale is not specified, it's based on the
    opacity at point. */
template <typename dust_model_type, typename grid_type>
mcrx::vec3d
mcrx::mlt_xfer<dust_model_type, grid_type>::
perturb_point(const vec3d& point, T_float scale) 
{
    // get 3 random numbers
  blitz::TinyVector<T_float, 3> r;
  this->rng().rnd(r);
  
  // see if scale was specified
  if(scale<0)
    scale = perturbation_scale(point, false);

  // draw radius from exp(-r/scale)/scale distribution. Note that this
  // implies that dP/dV = exp(-r/scale)/(4pi*scale*r^2) due to the
  // volume element, so we need to take that into account in the
  // probability. (Calculation double checked 110404.)
  const T_float radius = -log(r[0])*scale;

  DEBUG(1,std::cout << "perturbation scale " << scale << " size " << radius << std::endl;);
  // and the direction (from an isotropic distribution)
  const T_float dphi = r [1]*2*constants::pi; 
  const T_float dcos_theta = 2*r [2] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  const vec3d delta(radius*dsin_theta*cos(dphi),
		    radius*dsin_theta*sin(dphi),
		    radius*dcos_theta);
 
  const vec3d newpoint = point + delta;
  return newpoint;
}

/** Returns the perturbation scale based on the local opacity, or if
    use_image is set, also including perturbation in image plane.. */
template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
perturbation_scale(const vec3d& x, bool use_image) const
{
  const T_cell_tracker c = this->g.locate(x, false);
  T_float scale=0;
  if(c) {
    const T_densities rho (make_thread_local_copy(c.data()->get_absorber().densities()));
    if(any(rho>0)) {
	const T_float refkappa = 
	  this->dust_model.total_reference_abs_length(rho, this->b_);
      assert(refkappa>=0);
      scale = perturbation_scale_/refkappa;
      assert(scale>0);
    }
  }

  if(scale==0)
    scale = blitz::huge(T_float());

  if(use_image) {
    const T_float fov = sqrt((*this->current_emergence_->begin())->solid_angle());
    const vec3d campos = (*this->current_emergence_->begin())->position();
    const T_float imscale = image_perturbation_scale_* fov * mag(x-campos);
    DEBUG(1,std::cout << "Image perturbation, image scale " << imscale << " tau scale " << scale << std::endl;);
    scale = std::min(scale, imscale);
  }

  scale = std::min(scale, max_perturbation_scale_);
  assert(scale>0);
  return scale;
}

/** Returns probability (dP/dV) of perturbing x0 -> x, based on
    drawing radius from an exponential distribution and the angles
    from an isotropic one. See perturb_point for comments on the
    probability distribution. */
template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
perturbation_probability(const vec3d& x0, const vec3d& x, T_float scale) const
{
  assert(scale>0);
  const vec3d delta = x - x0;
    const T_float r=sqrt(dot(delta, delta));
  // see perturb_point for motivation of the p expression
  const T_float p=exp(-r/scale)/(4*constants::pi*scale*r*r);
  DEBUG(2,std::cout << r << '\t' << scale << '\t' << p << std::endl;);
  return p;
}
 
 template<typename T_iterator>
   void iter_add(T_iterator& i, int n) {assert(n>0);while(n--)++i;};
 template<typename T_iterator>
   void iter_sub(T_iterator& i, int n) {assert(n>0);while(n--)--i;};

/** Does a bidirectional mutation. Returns the perturbed path and the
    proposal ratio. */
template <typename dust_model_type, typename grid_type>
mcrx::transition_return_type<typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_lambda>
mcrx::mlt_xfer<dust_model_type, grid_type>::
bidir_mutation(const T_ray_path& p)
{
  DEBUG(1,std::cout << "\n*** Doing bidirectional mutation ***\n\n";);
  DEBUG(2,std::cout << "Starting path " << p << std::endl;);

  const int len=p.length();

  // Step 1: Draw a subpath to delete. This does not include endponts,
  // so can be between vertex 1 and n-2.

  // max # of vertices to delete (min is 0)
  const int nmax =len-4;

  // deletion weights
  const std::vector<T_float> dw = total_deletion_weight(p);
  // randomly draw a subpath
  const T_float dd=this->rng().rnd()*dw.back();
  // find which one it is by looping again
  int l,m;
  std::vector<T_float>::const_iterator i=dw.begin();
  for(l=1; l<=len-3; ++l)
    for(m=l+1; m<=len-2; ++m)
      if (*i++ > dd)
	goto hell;
 hell:
  // l,m should now be set to our deleted subpath
  DEBUG(1,std::cout << "Selected subpath " << l << "," << m << " for deletion\n";);
  T_float old_common; std::vector<T_float> p_edge;
  T_lambda old_i(path_intensity(p,old_common,p_edge,l,m));
  
  //number of vertices to delete
  const T_float nd = m-l-1;

  // Step 2: Draw a length of subpath to add.
  const T_float aa = 
    this->rng().rnd()*total_addition_weight(len, nd, 
					    (this->n_scat_min_==0));
  int na=0;
  T_float aw=addition_weight(len, nd, na, (this->n_scat_min_==0));
  while(aw<aa)
    aw += addition_weight(len, nd, ++na, (this->n_scat_min_==0));
  // na now contains the number of vertices to add. divide into ns and nt
  const int ns= int(floor(this->rng_.rnd()*(na+1)));
  assert(ns<=na);
  const int nc= na-ns;
  DEBUG(1,std::cout << "Selected " << ns << "," << nc << " vertices for addition\n";);

  // Step 4: Get the subpaths
  boost::shared_ptr<T_ray_path> s_subpath(new T_ray_path(p)); 
  T_ray_path c_subpath(s_subpath->split(l).split(m-l-2));
  //DEBUG(2,std::cout << "Source subpath: " << s_subpath << std::endl;);
  //DEBUG(2,std::cout << "Camera subpath: " << c_subpath << std::endl;);
  char buf[100];
  snprintf(buf,100,"BID %d %d %d %d",l,m,ns,nc);
  const std::string desc(buf);

  // Step 5: Do the actual mutation: extend these paths by ns and nc and merge.
  bool success = 
    extend_path(*s_subpath, ns, blitz::huge(int())) &&
    extend_path(c_subpath, nc, blitz::huge(int()));
  if(!success) {
    // extension attempt failed -- proposal failed
    DEBUG(1,std::cout << "Extension attempt failed\n";);
    return transition_return_type<T_lambda>(boost::shared_ptr<T_ray_path>(new T_ray_path()),
					    1,1,T_lambda(0.0*old_i), old_i, 1, desc);
  }
  merge_paths(*s_subpath, c_subpath);
  //DEBUG(2,std::cout << "New path: " << s_subpath << std::endl;);

  // Step 6: Evaluate result and probability ratio
  T_float new_common;
  T_lambda new_i(path_intensity(*s_subpath, new_common, p_edge, l, l+na+1));
  DEBUG(1,compact(compact(std::cout << "Old common=" << old_common << ", intensity: ", old_i) \
		  << "\nNew common=" << new_common << ", intensity: ", new_i) << std::endl;);

  // probability ratio is the ratio of deleting, adding, and extending
  T_float p_forward = 
    p_mutation(p)*
    deletion_weight(p, l, m)*
    addition_weight(len, nd, na, (this->n_scat_min_==0))*
    total_path_probability(*s_subpath, l, m-nd+na)/
    (dw.back()*total_addition_weight(len, nd, (this->n_scat_min_==0)));

  T_float p_backward = 
    p_mutation(*s_subpath)*
    deletion_weight(*s_subpath, l, m-nd+na)*
    addition_weight(len-nd+na, na, nd, (this->n_scat_min_==0))*
    total_path_probability(p, l, m) /
    (total_deletion_weight(*s_subpath).back()*
     total_addition_weight(len-nd+na, na, (this->n_scat_min_==0)));
  DEBUG(1,std::cout << "Forward p=" << p_forward << ", backward p=" << p_backward << std::endl;);
  assert(p_forward>0);
  assert(p_backward>=0 || sum(new_i)==0);

  // if we just deleted and added a single vertex, we need to add the
  // probability of perturbing that vertex, too
  if((nd==1)&&(na==1)) {
    assert(m-l==2);
    T_vlci vo = p.v_.begin();
    T_vlci vn = s_subpath->v_.begin();
    add_iter(vo, l+1);
    add_iter(vn, l+1);
    
    p_forward += 
      p_perturbation(p)*
      p_perturb_vertex(p,l+1)*
      perturbation_probability(vo->pos_, vn->pos_, 
			       perturbation_scale(vo->pos_, l==len-4));
    p_backward += 
      p_perturbation(*s_subpath)*
      p_perturb_vertex(*s_subpath,l+1)*
      perturbation_probability(vn->pos_, vo->pos_, 
			       perturbation_scale(vn->pos_, l==len-4));
  }
  if(!( p_forward>0 && (p_backward>=0 || sum(new_i)==0)))
    std::cout << *s_subpath << std::endl;
  assert(p_forward>0);
  assert(p_backward>=0 || sum(new_i)==0);

  return transition_return_type<T_lambda>(s_subpath, new_common, old_common,
					  new_i, old_i, 
					  p_backward/p_forward, desc);
}

  
/** Perturbs one of the scattering points. Returns the perturbed path
    and the proposal ratio. */
template <typename dust_model_type, typename grid_type>
mcrx::transition_return_type<typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_lambda>
mcrx::mlt_xfer<dust_model_type, grid_type>::
perturb_scattering_point(const T_ray_path& path)
{
  // number of perturbable points
  DEBUG(2,std::cout << "Perturbing a scattering point\n";);
  const int n = path.length()-4;
  
  assert(n>0);
  
  boost::shared_ptr<T_ray_path> new_path(new T_ray_path(path));
  
  // which scattering will be perturbed?
  
  // assign extra weight to perturbing the last point to improve image
  // distribution
  const T_float extra_weight = last_extra_weight_*n;
  int scattering = int(floor(this->rng().rnd()*(n+extra_weight)) + 2);

  if (scattering>n+1) {
    // n+1 is the last scattering vertex. if it points too far, it's
    // because of the extra weight
    scattering = n+1;
  }
  DEBUG(1,std::cout << "Perturbing vertex " << scattering <<  std::endl;);

  // set iterators pointing to the 3 vertices and 2 edges we are dealing with
  T_vli v2 = new_path->v_.begin();
  add_iter(v2, scattering-1);
  T_vli v0=v2++;
  T_vli v1=v2++;
  T_lli g2 = new_path->g_.begin();
  add_iter(g2, scattering-1);
  T_lli g1 = g2++;

  const T_float scale = perturbation_scale(v1->pos_, scattering==n+1);
  v1->pos_ = perturb_point(v1->pos_, scale);
  update_path_factors(*new_path, scattering-1, scattering+1);

  // we only need to evaluate the changed part of the path
  T_float old_common, new_common; std::vector<T_float> junk2;
  T_lambda old_i = path_intensity(path, old_common, junk2, 
				  scattering-1, scattering+1);
  T_lambda new_i = path_intensity(*new_path, new_common, junk2, 
				  scattering-1, scattering+1);
  assert(all(new_i==new_i));

  T_vlci ov1 = path.v_.begin();
  add_iter(ov1, scattering);
  const T_float new_scale = perturbation_scale(v1->pos_, scattering==n+1);

  const vec3d delta = v1->pos_ - ov1->pos_;
  DEBUG(1,std::cout << "Perturbing scattering point " << scattering << " from "
	<< ov1->pos_ << " to " << v1->pos_ << std::endl;);

  char buf[100];
  snprintf(buf,100,"MOV %d %d %5.3g %5.3g",
	   n, scattering, scale, sqrt(dot(delta,delta)));
  const std::string desc(buf);

  // if new_i==0, reject
  if((new_common==0) || (sum(new_i)==0))
    return transition_return_type<T_lambda>(new_path, new_common, old_common, 
					    new_i, old_i,
					    1.0, desc);

  // when calculating the probabilities, we need to also take into
  // account the probability of doing the same thing through a
  // bidirectional mutation
  const T_float p_forward = 
    // prob of perturbing path vertex to new_path vertex
    p_perturbation(path)*
    p_perturb_vertex(path,scattering)*
    perturbation_probability(ov1->pos_, v1->pos_, scale) +
    // prob of mutating path to new_path by deleting and adding a vertex
    p_mutation(path)*
    deletion_weight(path,scattering-1,scattering+1)*
    addition_weight(n+4, 1,1, (this->n_scat_min_==0))*
    total_path_probability(*new_path,scattering-1,scattering+1)/
    (total_deletion_weight(path).back()*
     total_addition_weight(n+4, 1, (this->n_scat_min_==0)));
  
  const T_float p_backward =
    // prob of perturbing new_path vertex to path vertex
    p_perturbation(*new_path)*
    p_perturb_vertex(*new_path,scattering)*
    perturbation_probability(v1->pos_, ov1->pos_, new_scale) +
    // prob of mutating new_path to path by deleting and adding a vertex
    p_mutation(*new_path)*
    deletion_weight(*new_path,scattering-1,scattering+1)*
    addition_weight(n+4, 1, 1, (this->n_scat_min_==0))*
    total_path_probability(path,scattering-1,scattering+1)/
    (total_deletion_weight(*new_path).back()*
     total_addition_weight(n+4, 1, (this->n_scat_min_==0)));

  if(!( p_forward>0 && p_backward>=0))
    std::cout << *new_path << std::endl;
  assert(p_forward>0);
  assert(p_backward>=0);

  return transition_return_type<T_lambda>(new_path, new_common, old_common, 
					  new_i, old_i,
					  p_backward/p_forward, desc);
}


/** Generates a proposal path and proposal ratio from an existing path. */
template <typename dust_model_type, typename grid_type>
mcrx::transition_return_type<typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_lambda>
mcrx::mlt_xfer<dust_model_type, grid_type>::
proposal(const T_ray_path& path)
{
  T_float p=this->rng().rnd()*total_weight(path);

  if (p<mutation_weight(path))
    return bidir_mutation(path);
  else
    return perturb_scattering_point(path);
}

/** Performs Metropolis sampling of the ray paths, starting with the
    initial path and adding contributions to the image. The
    probability of drawing the initial path p0 is used to weight all
    samples, which eliminates startup bias (according to the MLT
    paper). */
template <typename dust_model_type, typename grid_type>
int
mcrx::mlt_xfer<dust_model_type, grid_type>::
metropolis(const T_ray_path& start_path, int n, T_float w0, 
	   bool write_debug)
{
  
  std::ofstream debug;
  if(write_debug)
    debug.open("metro.debug");

  boost::shared_ptr<T_ray_path> current_path(new T_ray_path(start_path));
  T_float current_common; std::vector<T_float> p_edge;
  T_lambda current_i (make_thread_local_copy
		      (path_intensity(start_path, current_common, p_edge)));
  T_lambda temp_i(make_thread_local_copy(current_i));
  T_float temp_common;

  const T_float start_ii = sum(current_i);

  if(start_ii*current_common==0) {
    std::cout << "Initial path intensity=0, aborting this batch" << std::endl;
    return 0;
  }

  std::cout << "Initial path intensity: " << current_common*start_ii << std::endl;

  blitz::TinyVector<int,2> current_px =
    (*this->current_emergence_->begin())->project((++(++start_path.v_.rbegin()))->pos_);

  int rejections=0;
  bool underflow = any(current_i==0);

  for(int i=0; i<n; ++i) {
    // generate proposal
    DEBUG(1,std::cout << "\n\nGenerating proposal path" << std::endl;);
    assert(current_common>0);

    transition_return_type<T_lambda> prop = proposal(*current_path);    

    DEBUG(2,std::cout << "Proposal path is \n" << *prop.path << std::endl;);
    
    T_float i_ratio;
    if(!underflow) {
      // our sampled function is actually sum(i), so the ratio of new to
      // old is given by this formula
      assert(all(prop.i_old>0));
      assert(sum(current_i)>0);
      temp_i = current_i*prop.i_new/prop.i_old;
      // note that there's a missing current_common factor here...
      temp_common = (prop.new_common/prop.old_common);
      i_ratio = temp_common*sum(temp_i)/sum(current_i);
      // ... which we multiply in here
      temp_common *= current_common;
      assert(i_ratio>=0);
      assert(i_ratio<blitz::huge(T_float()));
    }
    else {
      // if there's underflow in the old intensity, we need to
      // evaluate the *entire* proposed path intensity, because we
      // can't take ratios
      temp_i = path_intensity(*prop.path, temp_common, p_edge);
      i_ratio = (temp_common/current_common)*sum(temp_i)/sum(current_i);
      assert(i_ratio>=0);
      assert(i_ratio<blitz::huge(T_float()));
    }
     
    // project final scattering point on image (assume only one camera)
    blitz::TinyVector<int,2> px = 
      (*this->current_emergence_->begin())->project((++(++prop.path->v_.rbegin()))->pos_);
    
    if(px[0]==-1)
      i_ratio=0;
  

    // Generate acceptance probability. There's an inherent problem
    // with calculating probabilities when the proposal intensity is
    // zero, so in that case we just set a to 0.
    const T_float a = (i_ratio>0) ? std::min(prop.p_ratio * i_ratio , 1.0) : 0;
    
    assert(a==a);

    DEBUG(1,std::cout << "Proposal path: " << prop.path << " int ratio=" \
	  << i_ratio << " trans ratio=" << prop.p_ratio \
	  << " a=" << a << std::endl;);

    // add the conditional path to the image. What we add only depends
    // on the ratios of wavelengths, so is independent of the common
    // factors
    if( (a>0) &&
	(prop.path->length()-4 >= this->n_scat_min_) &&
	(prop.path->length()-4 <= this->n_scat_max_) )
      (*this->current_emergence_->begin())->add(px, a*w0, temp_i/sum(temp_i));
    
    if( (a<1) &&
	(current_path->length()-4 >= this->n_scat_min_) &&
	(current_path->length()-4 <= this->n_scat_max_) )
      (*this->current_emergence_->begin())->add(current_px, (1-a)*w0,
						current_i/sum(current_i));
    
    std::string desc;

    if ( (a>=1) || (this->rng().rnd()<a) ) {
      // we accept
      current_path = prop.path;
      if(i%100==0 && !underflow) {
	// regenerate full path so we don't get weird numeric issues 
	current_i=path_intensity(*current_path, current_common, p_edge);
      }
      else {
	current_i = temp_i;
	current_common = temp_common;
      }

      current_px = px;
      DEBUG(1,compact(std::cout << "Proposal accepted, current common=" << current_common \
		      << ", i=", current_i) << " @ " << current_px << std::endl;);
      desc = "PROP A ";
      // it is POSSIBLE to have zero intensity here due to underflow
      // with an accepted proposal that lowered intensity. we need to
      // deal with that the next time around
      if(underflow && all(current_i>0))
	std::cout << "Recovered from intensity underflow" << std::endl;
      underflow = any(current_i==0);

      if(underflow) {
	std::cout << "Warning: intensity underflow for accepted path in mutation " << i << '\n';
      DEBUG(1, compact(cout << "Current intensity: ", current_i) << "\n\nExisting path is \n" << *current_path << std::endl;);
      }
    }
    else {
      DEBUG(1,compact(std::cout << "Proposal rejected, current i=", current_i) \
	    << " @ " << current_px << std::endl;);
      desc = "PROP R ";
      ++rejections;
    }

    if(write_debug) {
      // we do this because lexical_cast doesn't support formatting...
      char buf[100];
      snprintf(buf,100,"%9.3e %9.3e %9.3e ",
	       i_ratio, prop.p_ratio, sum(current_i));
      desc=desc+buf+prop.desc;
      debug << desc << '\n';
    }
  }

  std::cout << "Acceptance ratio " << 1.0-1.0*rejections/n << std::endl;

  return rejections;
}

template <typename dust_model_type, typename grid_type>
typename mcrx::mlt_xfer<dust_model_type, grid_type>::T_ray_path
mcrx::mlt_xfer<dust_model_type, grid_type>::
make_path(const std::vector<vec3d>& pos)
{
  T_ray_path p(T_ray_path::create_empty(true, 3));
  for(int i=0;i<pos.size()-1;++i) {
    T_vertex v(3);
    v.pos_=pos[i];
    v.s_=1;v.e_=1;
    p.v_.push_back(v);
    if(i<pos.size()-2)
      p.g_.push_back(T_lambda(3));
  }
  T_ray_path pp(T_ray_path::create_empty(false, 3));
  T_vertex v(3);
  v.pos_=pos.back();
  v.s_=1;v.e_=1;
  pp.v_.push_front(v);
  //std::cout << p << pp << std::endl;
  merge_paths(p,pp);
  update_path_factors(p,1,p.length()-2);
  return p;
}

template <typename dust_model_type, typename grid_type>
void
mcrx::mlt_xfer<dust_model_type, grid_type>::
pathtest()
{
  std::vector<vec3d> v;
  v.push_back(vec3d(0.01,0.01,-0.01));
  //v.push_back(vec3d(-1.15823,-0.803942,-0.165162));
  v.push_back(vec3d(-0.731713,-1.57416,-0.184028));
  v.push_back(vec3d(-0.697923,-1.6258,-0.0484195));
  v.push_back(vec3d(46172.4,19141.5,1308.85));
  T_ray_path p(make_path(v));
  std::cout << p << '\n';
  std::cout << path_intensity(p) << path_probability(p,1,p.length()-3) << std::endl;
}

template <typename dust_model_type, typename grid_type>
void
mcrx::mlt_xfer<dust_model_type, grid_type>::
bidirectional_shoot(int ns_max, int nc_max)
{
  generate_bidir_paths(true, (this->n_scat_min_>0), ns_max, nc_max);
}
/*
  //  pathtest();return 0;
  int n=0;
  for(int i=0; i<nrays; ++i) {
    n+= 
  }
  std::cout << nrays << " bidirectional path pairs generated " << n << " samples"<<std::endl;

  return nrays;
}
*/

/** Runs n1 iterations, where it generates a random ray path based on
    the forward probabilities and then uses this as a seed for n2
    metropolis trials. 

    The number of scatterings in the seed path is determined randomly
    before calling shoot, because the peel-off technique means we
    don't really have a set number of scatterings. (We could use when
    it exits the volume, but then the last scattering point might not
    connect in any meaningful way to the observer anyway.) No, for now
    we use the exit number of scatterings.  */
template <typename dust_model_type, typename grid_type>
mcrx::T_float
mcrx::mlt_xfer<dust_model_type, grid_type>::
metropolis_driver(int nseed, int nmetropolis, int nmutations)
{
  max_n_track_=-1;

  std::cout << "Generating " << nseed << " seed paths" << std::endl;

  // generate seed paths
  std::map<int, typename T_rng::T_state> seed_states;
  for(int i=0; i<nseed; ++i) {
    // save the state so we can recover
    seed_states[seed_weights_.size()]=this->rng().get_generator().getState();

    generate_bidir_paths(false, (this->n_scat_min_>0) );
  }
  // save final state
  typename T_rng::T_state final= this->rng().get_generator().getState();

  const int n=seed_weights_.size();
  std::cout << nseed << " seed shoots generated " << n << " seed paths:\n";

  assert(seed_lengths_.size()==n);

  /*
  //now go over paths and calcualte weights and intensities
  std::vector<T_float> seed_weights;
  T_float junk1; vector<T_float> junk2;
  for(int i=0; i<n; ++i) {    
    // for very long paths, the probability can underflow to zero. In
    // that case we can't define the weight, so we just ignore it by
    // setting weight to 0
    T_lambda cur_i(path_intensity(seed_paths_[i], junk1, junk2));
    // note that we divide by n_lambda here to get the normalization right
    T_float cur_f = any(cur_i==0) ? 0.0 : sum(cur_i);
    T_float cur_weight = (seed_probs_[i]>0) ? (cur_f/seed_probs_[i]) : 0.0;
    seed_weights.push_back(cur_weight + 
			   ((i==0) ? 0.0 : seed_weights.back()));
    //printf("Seed path %d with p=%3g, i=%3g, w=%3g\n", i, seed_probs_[i], cur_f, cur_weight);
    //std::cout << seed_paths_[i] << std::endl;
    assert(cur_f==0 || all(cur_i>0));
    assert(cur_f>=0);
    assert(cur_weight<blitz::huge(T_float()));    
  }
  */

  // now we resample the seed paths and use those for MLT runs with
  // the average weight as the initial weight. Note that it SHOULD be
  // nseed here.
  const T_float mean_weight=seed_weights_.back()/nseed;
  std::cout << "Mean path weight was " << mean_weight << std::endl;

  // Note that to get the intensity we would have gotten for the seed
  // paths, we have to divide by nSEED, not n, because the paths were
  // split up for each successive scattering here while they count as
  // one ray in normal runs.

  int rej=0;
  for(int i=0; i<nmetropolis; ++i) {
    // pick path p randomly in cumulative distro
    const int p =
      std::lower_bound (seed_weights_.begin(), seed_weights_.end(), 
			seed_weights_.back()*this->rng().rnd()) - //i/(nseed-1)) -
      seed_weights_.begin();
    assert(p>=0);
    assert(p<n);

    // now we need to regenerate the path. first find the state
    typename std::map<int, typename T_rng::T_state>::const_iterator st =
      seed_states.lower_bound (p);
    if(st->first>p) --st;
    assert(st->first<=p);
    final=this->rng().get_generator().getState();
    this->rng().get_generator().setState(st->second);

    // now recreate the path
    T_ray_path current_path(recreate_bidir_path(seed_lengths_[p]));
    // reset the generator
    this->rng().get_generator().setState(final);
    
    std::cout << "\n\nPicked seed path " << p <<" with weight " 
	 << seed_weights_[p]-(p==0 ? 0.0 : seed_weights_[p-1]) << ":\n" 
	 << current_path << std::endl;

    // since we use the same weight, we can multiply by the mean weight later
    rej+= metropolis(current_path, nmutations, 1, false);
  }
  T_float ar = 1.0-1.0*rej/(nmutations*nmetropolis);
  std::cout << "Overall acceptance ratio " << ar << std::endl;

  return mean_weight;
}


/** Constructor. Sets up the transfer class to handle a specified grid
    and send the output to a specified emergence object. */
template <typename dust_model_type, typename grid_type>
mcrx::mlt_xfer<dust_model_type, grid_type>::
mlt_xfer(T_grid& grid, const T_emission& emi, T_emergence& eme, 
	 const T_dust_model& dm, T_rng_policy& rng,
	 const T_biaser& b, int minscat, int maxscat) :
  xfer<dust_model_type, grid_type>
  (grid, emi, eme, dm, rng, b, 1e-4, 1e2, false, false, minscat, maxscat)
{
  perturbation_scale_=0.5;
  max_perturbation_scale_=100.0;
  mutation_weight_=0.5;
  perturbation_weight_=0.5;
  image_perturbation_scale_=0.1;
  last_extra_weight_ = 0.4;

  // get distributions (XXX these should be in the path but for now this'll do)
  T_float junk;
  ray_distribution_pair<T_lambda> e = 
    this->emission_.emit(this->b_, this->rng_, junk);
  emdistr_ =  e.distribution();
  camdistr_ = (*this->current_emergence_->begin())->distribution();
  dustdistr_ = boost::shared_ptr<angular_distribution<T_lambda> >
    (this->dust_model.get_scatterer(0).phase_function().clone());
  // force camera to allocate 
  for(typename T_emergence::iterator i=this->current_emergence_->begin();
      i!=this->current_emergence_->end(); ++i)
    (*i)->allocate(this->emission_.zero_lambda());

  consolidate_arrays(camera_intensity_);
}


//!  Copy constructor
/*!  Makes a copy of an xfer object, including the current cell, lost
  energy and random number generator state.*/
template <typename dust_model_type, typename grid_type>
mcrx::mlt_xfer<dust_model_type, grid_type>::
mlt_xfer(const mlt_xfer& rhs):
  xfer<dust_model_type, grid_type>(rhs),
  perturbation_scale_(rhs.perturbation_scale_),
  max_perturbation_scale_(rhs.max_perturbation_scale_),
  mutation_weight_(rhs.mutation_weight_),
  perturbation_weight_(rhs.perturbation_weight_),
  image_perturbation_scale_(rhs.image_perturbation_scale_),
  last_extra_weight_(rhs.last_extra_weight_),
  max_n_track_(rhs.max_n_track_)
{
  // get distributions (XXX these should be in the path but for now this'll do)
  T_float junk;
  ray_distribution_pair<T_lambda> e = 
    this->emission_.emit(this->b_, this->rng(), junk);
  emdistr_ =  e.distribution();
  camdistr_ = (*this->current_emergence_->begin())->distribution();
  dustdistr_ = boost::shared_ptr<angular_distribution<T_lambda> >(this->dust_model.get_scatterer(0).phase_function().clone());
  // force camera to allocate 
  for(typename T_emergence::iterator i=this->current_emergence_->begin();
      i!=this->current_emergence_->end(); ++i)
    (*i)->allocate(this->emission_.zero_lambda());

  consolidate_arrays(camera_intensity_);
}

/** We need to redo this in the mlt_xfer class to add the
    camera_intensity_ array. */
template <typename dust_model_type, typename grid_type>
void mcrx::mlt_xfer<dust_model_type, grid_type>::
consolidate_arrays(array_1&)
{
  const int narr=11;
  /// Pad the allocated memory so different threads don't share cache
  /// lines
  const int buffer=2; 
  array_2 temp(narr+buffer*2, this->emission_.zero_lambda().size());
  threadLocal_warn(temp);
  int i=buffer;
  /// Because temp is local, we can't return *weak* references.
  this->norm_.reference(temp(i++,blitz::Range::all()));
  this->i_in_.reference(temp(i++,blitz::Range::all()));
  this->i_out_.reference(temp(i++,blitz::Range::all()));
  this->dtau_.reference(temp(i++,blitz::Range::all()));
  this->tau_exit_.reference(temp(i++,blitz::Range::all()));
  this->onemexptau_exit_.reference(temp(i++,blitz::Range::all()));
  this->scattering_tau_.reference(temp(i++,blitz::Range::all()));
  this->original_intensity_.reference(temp(i++,blitz::Range::all()));
  this->intensity_to_camera_.reference(temp(i++,blitz::Range::all()));
  camera_intensity_.reference(temp(i++,blitz::Range::all()));
  assert(i==narr+buffer-1);
}

#endif



